## 
## The downloaded binary packages are in
##  /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpeU3KYM/downloaded_packages
## 
## The downloaded binary packages are in
##  /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpeU3KYM/downloaded_packages
## 
## The downloaded binary packages are in
##  /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpeU3KYM/downloaded_packages
## 
## The downloaded binary packages are in
##  /var/folders/7n/x74qctp91rng390gx0z9hmd80000gn/T//RtmpeU3KYM/downloaded_packages
load(here("jk_code", "JK_remove_mtrp.rds"))


 SO4@meta.data <- SO4@meta.data %>% 
  mutate(subclass_MD = dplyr::case_when(
    seurat_clusters == 0  ~ "type_1",
    seurat_clusters == 1  ~ "type_1",
    seurat_clusters == 2  ~ "type_3",
    seurat_clusters == 3  ~ "type_2",
    seurat_clusters == 4  ~ "type_4",
    seurat_clusters == 5  ~ "type_5",

  ))

SO4@meta.data$subclass_MD <- factor(SO4@meta.data$subclass_MD , levels = c("type_1", "type_2", "type_3", "type_4","type_5"))

Idents(SO4) <- SO4@meta.data$subclass_MD

DimPlot(object = SO4, reduction = "umap", group.by = "subclass_MD", label = TRUE)

DimPlot(object = SO4, reduction = "umap", label = TRUE)

Idents(SO4) <- "subclass_MD"

DimPlot(SO4,split.by = "sample",group.by = "seurat_clusters")

# Find markers for all clusters (adjust parameters as needed)

Allmarkers <- FindAllMarkers(SO4, only.pos = TRUE)
## Calculating cluster type_1
## Warning: The `slot` argument of `GetAssayData()` is deprecated as of SeuratObject 5.0.0.
## ℹ Please use the `layer` argument instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## Warning: `PackageCheck()` was deprecated in SeuratObject 5.0.0.
## ℹ Please use `rlang::check_installed()` instead.
## ℹ The deprecated feature was likely used in the Seurat package.
##   Please report the issue at <https://github.com/satijalab/seurat/issues>.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.
## For a (much!) faster implementation of the Wilcoxon Rank Sum Test,
## (default method for FindMarkers) please install the presto package
## --------------------------------------------
## install.packages('devtools')
## devtools::install_github('immunogenomics/presto')
## --------------------------------------------
## After installation of presto, Seurat will automatically use the more 
## efficient implementation (no further action necessary).
## This message will be shown once per session
## Calculating cluster type_2
## Calculating cluster type_3
## Calculating cluster type_4
## Calculating cluster type_5
Allmarkers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1)
Allmarkers %>%
    group_by(cluster) %>%
    dplyr::filter(avg_log2FC > 1) %>%
    slice_head(n = 10) %>%
    ungroup() -> top10
DoHeatmap(SO4, features = top10$gene) + NoLegend()
## Warning in DoHeatmap(SO4, features = top10$gene): The following features were
## omitted as they were not found in the scale.data slot for the SCT assay: Phkg1,
## Gm26982, Rtp4, Rsad2, Oasl2, Oasl1, Gstm6

type_1_markers <- Allmarkers[Allmarkers$cluster == "type_1", ]
type_2_markers <- Allmarkers[Allmarkers$cluster == "type_2", ]
type_3_markers <- Allmarkers[Allmarkers$cluster == "type_3", ]
type_4_markers <- Allmarkers[Allmarkers$cluster == "type_4", ]
type_5_markers <- Allmarkers[Allmarkers$cluster == "type_5", ]
markers.to.plot1 <- c("Ramp3",
                      "Slc12a1",
                      "Egf",      # PT-S1
                      "Fos",
                      "Cxcl10")      # PT-S3
                
DotPlot(SO4,
features = markers.to.plot1,
dot.scale = 8,
dot.min = 0,
scale.max = 100,
scale.min = 0,
col.min = -2.5,
col.max = 2.5)+
coord_flip()

1 Type 1

df<- type_1_markers %>% arrange(desc(avg_log2FC))

df2 <- df %>% filter(p_val_adj < 0.05)

DEG_list <- df2

markers <- DEG_list %>% rownames_to_column(var="SYMBOL")

head(markers, n = 50)
ENTREZ_list <- bitr(
  geneID = rownames(DEG_list),
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Mm.eg.db
)
## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 3.52% of input gene IDs are fail to map...
markers <-  ENTREZ_list %>% inner_join(markers, by = "SYMBOL")

# Removing genes that are not statistically significant. 
markers <-  markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)

pos.markers <-  markers %>% dplyr::filter(avg_log2FC > 0) %>%  arrange(desc(abs(avg_log2FC)))
#head(pos.markers, n = 50)

pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)

pos_go <- enrichGO(gene = pos.ranks,           #a vector of entrez gene id
                   OrgDb = "org.Mm.eg.db",    
                   ont = "BP",
                   readable = TRUE)              #whether mapping gene ID to gene Name

pos_go
## #
## # over-representation test
## #
## #...@organism     Mus musculus 
## #...@ontology     BP 
## #...@keytype      ENTREZID 
## #...@gene     chr [1:2714] "102633582" "69480" "102632087" "11551" "73955" "14062" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...1077 enriched terms found
## 'data.frame':    1077 obs. of  12 variables:
##  $ ID            : chr  "GO:0008380" "GO:0006457" "GO:0034976" "GO:0048193" ...
##  $ Description   : chr  "RNA splicing" "protein folding" "response to endoplasmic reticulum stress" "Golgi vesicle transport" ...
##  $ GeneRatio     : chr  "127/2524" "74/2524" "88/2524" "86/2524" ...
##  $ BgRatio       : chr  "464/28928" "185/28928" "279/28928" "293/28928" ...
##  $ RichFactor    : num  0.274 0.4 0.315 0.294 0.283 ...
##  $ FoldEnrichment: num  3.14 4.58 3.61 3.36 3.24 ...
##  $ zScore        : num  14.3 15.1 13.6 12.6 12.3 ...
##  $ pvalue        : num  2.57e-32 6.54e-31 1.49e-27 1.56e-24 4.66e-24 ...
##  $ p.adjust      : num  1.57e-28 1.99e-27 3.02e-24 2.37e-21 4.06e-21 ...
##  $ qvalue        : num  1.17e-28 1.49e-27 2.26e-24 1.77e-21 3.04e-21 ...
##  $ geneID        : chr  "Zfp64/Khdrbs2/Prmt1/Ess2/Sf3b4/Snrpn/Cirbp/Thoc3/Dcps/Thoc6/Rbm42/Ppil1/U2af1l4/Rbmx2/Syf2/Ppih/Snrpa/BC005624/"| __truncated__ "Prdx4/Hspb6/H2-DMb1/Ppic/Sdf2l1/Elp6/Ric3/Clu/B2m/Qsox2/Ppil1/Pdia4/Ppih/Sdf2/Qsox1/Arl2/Fkbp9/Aip/Tor1a/Nudcd2"| __truncated__ "Atp2a3/Nupr1/Fbxo2/Nr1h3/Rnf186/Serp2/Sdf2l1/Fbxo6/Flot1/Marcks/Atg10/Clu/Pdia4/Creb3/Atf6b/Pdia6/Tor1a/H13/Stu"| __truncated__ "Tmed6/Atp9a/Stx18/Lyplal1/Gosr2/Krt18/Tmed4/Ergic3/Rab13/Cog8/Tmed1/Tmed3/Trappc4/Ccdc22/Golga7/Tmed10/Tmed9/Go"| __truncated__ ...
##  $ Count         : int  127 74 88 86 89 89 89 107 90 46 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
 dotplot(pos_go) +
    ggtitle("") +
    theme_classic() + 
    theme(
        plot.title = element_text(hjust = 0.5),
        legend.position = "left",
        axis.text.y = element_text(hjust = 0, size = 10)) +
    scale_y_discrete(position = "right", 
                     labels = function(x) str_wrap(x, width = 25))  # Wrap y-axis labels to 2 lines
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

# Type 2

  df<- type_2_markers %>% arrange(desc(avg_log2FC))


df2 <- df %>% filter(p_val_adj < 0.05)

DEG_list <- df2

markers <- DEG_list %>% rownames_to_column(var="SYMBOL")

head(markers, n = 50)
ENTREZ_list <- bitr(
  geneID = rownames(DEG_list),
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Mm.eg.db
)
## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 22.64% of input gene IDs are fail to map...
markers <-  ENTREZ_list %>% inner_join(markers, by = "SYMBOL")

# Removing genes that are not statistically significant. 
markers <-  markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)

pos.markers <-  markers %>% dplyr::filter(avg_log2FC > 0) %>%  arrange(desc(abs(avg_log2FC)))
#head(pos.markers, n = 50)

pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)

pos_go <- enrichGO(gene = pos.ranks,           #a vector of entrez gene id
                   OrgDb = "org.Mm.eg.db",    
                   ont = "BP",
                   readable = TRUE)              #whether mapping gene ID to gene Name

pos_go
## #
## # over-representation test
## #
## #...@organism     Mus musculus 
## #...@ontology     BP 
## #...@keytype      ENTREZID 
## #...@gene     chr [1:622] "100503991" "100416292" "73010" "102638255" "791329" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...828 enriched terms found
## 'data.frame':    828 obs. of  12 variables:
##  $ ID            : chr  "GO:0008380" "GO:1903311" "GO:0043484" "GO:0000375" ...
##  $ Description   : chr  "RNA splicing" "regulation of mRNA metabolic process" "regulation of RNA splicing" "RNA splicing, via transesterification reactions" ...
##  $ GeneRatio     : chr  "62/594" "46/594" "33/594" "40/594" ...
##  $ BgRatio       : chr  "464/28928" "320/28928" "194/28928" "315/28928" ...
##  $ RichFactor    : num  0.134 0.144 0.17 0.127 0.127 ...
##  $ FoldEnrichment: num  6.51 7 8.28 6.18 6.18 ...
##  $ zScore        : num  17.3 15.6 14.7 13.4 13.4 ...
##  $ pvalue        : num  5.33e-32 2.48e-25 7.92e-21 3.65e-20 3.65e-20 ...
##  $ p.adjust      : num  2.29e-28 5.33e-22 1.14e-17 2.61e-17 2.61e-17 ...
##  $ qvalue        : num  1.60e-28 3.72e-22 7.94e-18 1.83e-17 1.83e-17 ...
##  $ geneID        : chr  "Paxbp1/Malat1/Ddx17/Supt6/Setx/Hnrnph1/Prpf38b/Rbm12b2/Sfpq/Zfp638/Rnpc3/Atxn7/Prpf4b/Prpf39/Rbm6/Clk4/Hnrnpa2b"| __truncated__ "Malat1/Ddx17/Ago3/Supt6/Pan3/Trim71/Rbm33/Tut4/Hnrnpa2b1/Cnot6l/Tardbp/Ago2/Rc3h2/Tnrc6a/Tnrc6b/Rbm15/Hnrnpu/Ac"| __truncated__ "Malat1/Ddx17/Setx/Hnrnph1/Rbm12b2/Zfp638/Atxn7/Clk4/Hnrnpa2b1/Clk1/Cdk12/Rbm15/Hnrnpu/Acin1/Ythdc1/Tra2a/Rbm5/S"| __truncated__ "Paxbp1/Malat1/Ddx17/Setx/Sfpq/Rnpc3/Prpf4b/Prpf39/Rbm6/Hnrnpa2b1/Srrm2/Luc7l2/Cdk13/Rbm15/Hnrnpu/Acin1/Ythdc1/T"| __truncated__ ...
##  $ Count         : int  62 46 33 40 40 40 29 34 29 29 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
 dotplot(pos_go) +
    ggtitle("") +
    theme_classic() + 
    theme(
        plot.title = element_text(hjust = 0.5),
        legend.position = "left",
        axis.text.y = element_text(hjust = 0, size = 10)) +
    scale_y_discrete(position = "right", 
                     labels = function(x) str_wrap(x, width = 25))  # Wrap y-axis labels to 2 li
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

2 Type 3

df<- type_3_markers %>% arrange(desc(avg_log2FC))

df
df2 <- df %>% filter(p_val_adj < 0.05)

DEG_list <- df2

markers <- DEG_list %>% rownames_to_column(var="SYMBOL")

head(markers, n = 50)
ENTREZ_list <- bitr(
  geneID = rownames(DEG_list),
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Mm.eg.db
)
## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 14.19% of input gene IDs are fail to map...
markers <-  ENTREZ_list %>% inner_join(markers, by = "SYMBOL")

# Removing genes that are not statistically significant. 
markers <-  markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)

pos.markers <-  markers %>% dplyr::filter(avg_log2FC > 0) %>%  arrange(desc(abs(avg_log2FC)))
#head(pos.markers, n = 50)

pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)

pos_go <- enrichGO(gene = pos.ranks,           #a vector of entrez gene id
                   OrgDb = "org.Mm.eg.db",    
                   ont = "BP",
                   readable = TRUE)              #whether mapping gene ID to gene Name

pos_go
## #
## # over-representation test
## #
## #...@organism     Mus musculus 
## #...@ontology     BP 
## #...@keytype      ENTREZID 
## #...@gene     chr [1:1179] "213417" "330428" "16513" "268379" "53315" "15220" "69824" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...940 enriched terms found
## 'data.frame':    940 obs. of  12 variables:
##  $ ID            : chr  "GO:0015980" "GO:0045333" "GO:0009060" "GO:0006119" ...
##  $ Description   : chr  "energy derivation by oxidation of organic compounds" "cellular respiration" "aerobic respiration" "oxidative phosphorylation" ...
##  $ GeneRatio     : chr  "131/1086" "114/1086" "99/1086" "75/1086" ...
##  $ BgRatio       : chr  "380/28928" "271/28928" "206/28928" "154/28928" ...
##  $ RichFactor    : num  0.345 0.421 0.481 0.487 0.512 ...
##  $ FoldEnrichment: num  9.18 11.21 12.8 12.97 13.63 ...
##  $ zScore        : num  31.7 33.3 33.6 29.4 28.2 ...
##  $ pvalue        : num  1.38e-89 2.39e-89 1.08e-84 8.49e-65 4.31e-58 ...
##  $ p.adjust      : num  6.21e-86 6.21e-86 1.87e-81 1.11e-61 4.49e-55 ...
##  $ qvalue        : num  4.59e-86 4.59e-86 1.38e-81 8.17e-62 3.32e-55 ...
##  $ geneID        : chr  "Mt3/Adrb1/Mrap2/Ppargc1a/Gabarapl1/Idh2/Cycs/Idh1/Cox5b/Got1/Eva1a/Cox7b/Idh3a/Ppp1r3b/Cox6c/Uqcrq/Chchd10/Uqcr"| __truncated__ "Ppargc1a/Idh2/Cycs/Idh1/Cox5b/Got1/Cox7b/Idh3a/Cox6c/Uqcrq/Chchd10/Uqcr10/Aco2/Cox5a/Uqcr11/Ndufb8/Slc25a12/Cox"| __truncated__ "Idh2/Cycs/Idh1/Cox5b/Cox7b/Idh3a/Cox6c/Uqcrq/Chchd10/Uqcr10/Aco2/Cox5a/Uqcr11/Ndufb8/Cox6a1/Cox6b1/Cs/Ndufs2/Pp"| __truncated__ "Cycs/Cox5b/Cox7b/Cox6c/Uqcrq/Chchd10/Uqcr10/Cox5a/Uqcr11/Ndufb8/Cox6a1/Cox6b1/Ndufs2/Ppif/Cyc1/Uqcrfs1/Cox7a2/C"| __truncated__ ...
##  $ Count         : int  131 114 99 75 65 63 56 48 48 52 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
 dotplot(pos_go) +
    ggtitle("") +
    theme_classic() + 
    theme(
        plot.title = element_text(hjust = 0.5),
        legend.position = "left",
        axis.text.y = element_text(hjust = 0, size = 10)) +
    scale_y_discrete(position = "right", 
                     labels = function(x) str_wrap(x, width = 25))  # Wrap y-axis labels to 2 lin
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

3 Type 4

df<- type_4_markers %>% arrange(desc(avg_log2FC))

df
df2 <- df %>% filter(p_val_adj < 0.05)

DEG_list <- df2

markers <- DEG_list %>% rownames_to_column(var="SYMBOL")

head(markers, n = 50)
ENTREZ_list <- bitr(
  geneID = rownames(DEG_list),
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Mm.eg.db
)
## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 37.44% of input gene IDs are fail to map...
markers <-  ENTREZ_list %>% inner_join(markers, by = "SYMBOL")

# Removing genes that are not statistically significant. 
markers <-  markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)

pos.markers <-  markers %>% dplyr::filter(avg_log2FC > 0) %>%  arrange(desc(abs(avg_log2FC)))
#head(pos.markers, n = 50)

pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)

pos_go <- enrichGO(gene = pos.ranks,           #a vector of entrez gene id
                   OrgDb = "org.Mm.eg.db",    
                   ont = "BP",
                   readable = TRUE)              #whether mapping gene ID to gene Name

pos_go
## #
## # over-representation test
## #
## #...@organism     Mus musculus 
## #...@ontology     BP 
## #...@keytype      ENTREZID 
## #...@gene     chr [1:122] "16149" "13654" "23833" "214855" "14282" "14825" "13653" ...
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...741 enriched terms found
## 'data.frame':    741 obs. of  12 variables:
##  $ ID            : chr  "GO:0010563" "GO:0045936" "GO:0042326" "GO:0001933" ...
##  $ Description   : chr  "negative regulation of phosphorus metabolic process" "negative regulation of phosphate metabolic process" "negative regulation of phosphorylation" "negative regulation of protein phosphorylation" ...
##  $ GeneRatio     : chr  "18/101" "18/101" "17/101" "16/101" ...
##  $ BgRatio       : chr  "388/28928" "388/28928" "337/28928" "316/28928" ...
##  $ RichFactor    : num  0.0464 0.0464 0.0504 0.0506 0.0393 ...
##  $ FoldEnrichment: num  13.3 13.3 14.4 14.5 11.3 ...
##  $ zScore        : num  14.4 14.4 14.7 14.3 13.1 ...
##  $ pvalue        : num  1.80e-15 1.80e-15 2.96e-15 1.95e-14 3.13e-14 ...
##  $ p.adjust      : num  2.27e-12 2.27e-12 2.48e-12 1.23e-11 1.58e-11 ...
##  $ qvalue        : num  1.30e-12 1.30e-12 1.42e-12 7.04e-12 9.06e-12 ...
##  $ geneID        : chr  "Socs3/Dusp1/Jun/Gadd45b/Ppp1r15a/Irf1/Hspb1/Socs1/Errfi1/Cdkn1a/Trib1/Zc3h12a/Ddit4/Gadd45a/Irs2/Dnaja1/Taf7/Actb" "Socs3/Dusp1/Jun/Gadd45b/Ppp1r15a/Irf1/Hspb1/Socs1/Errfi1/Cdkn1a/Trib1/Zc3h12a/Ddit4/Gadd45a/Irs2/Dnaja1/Taf7/Actb" "Socs3/Dusp1/Jun/Gadd45b/Ppp1r15a/Irf1/Hspb1/Socs1/Errfi1/Cdkn1a/Trib1/Zc3h12a/Gadd45a/Irs2/Dnaja1/Taf7/Actb" "Socs3/Dusp1/Jun/Gadd45b/Ppp1r15a/Irf1/Hspb1/Socs1/Errfi1/Cdkn1a/Trib1/Zc3h12a/Gadd45a/Dnaja1/Taf7/Actb" ...
##  $ Count         : int  18 18 17 16 18 17 7 12 7 14 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
 dotplot(pos_go) +
    ggtitle("") +
    theme_classic() + 
    theme(
        plot.title = element_text(hjust = 0.5),
        legend.position = "left",
        axis.text.y = element_text(hjust = 0, size = 10)) +
    scale_y_discrete(position = "right", 
                     labels = function(x) str_wrap(x, width = 25))  # Wrap y-axis labels to 2 li
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.

4 Type 5

df<- type_5_markers %>% arrange(desc(avg_log2FC))

df
df2 <- df %>% filter(p_val_adj < 0.05)

DEG_list <- df2

markers <- DEG_list %>% rownames_to_column(var="SYMBOL")

head(markers, n = 50)
ENTREZ_list <- bitr(
  geneID = rownames(DEG_list),
  fromType = "SYMBOL",
  toType = "ENTREZID",
  OrgDb = org.Mm.eg.db
)
## 'select()' returned 1:1 mapping between keys and columns
## Warning in bitr(geneID = rownames(DEG_list), fromType = "SYMBOL", toType =
## "ENTREZID", : 71.43% of input gene IDs are fail to map...
markers <-  ENTREZ_list %>% inner_join(markers, by = "SYMBOL")

# Removing genes that are not statistically significant. 
markers <-  markers %>% dplyr::filter(p_val_adj < 0.05)
#head(markers, n = 50)

pos.markers <-  markers %>% dplyr::filter(avg_log2FC > 0) %>%  arrange(desc(abs(avg_log2FC)))
#head(pos.markers, n = 50)

pos.ranks <- pos.markers$ENTREZID[abs(pos.markers$avg_log2FC) > 0]
#head(pos.ranks)

pos_go <- enrichGO(gene = pos.ranks,           #a vector of entrez gene id
                   OrgDb = "org.Mm.eg.db",    
                   ont = "BP",
                   readable = TRUE)              #whether mapping gene ID to gene Name

pos_go
## #
## # over-representation test
## #
## #...@organism     Mus musculus 
## #...@ontology     BP 
## #...@keytype      ENTREZID 
## #...@gene     chr [1:6] "231655" "18682" "58185" "15957" "23962" "67775"
## #...pvalues adjusted by 'BH' with cutoff <0.05 
## #...79 enriched terms found
## 'data.frame':    79 obs. of  12 variables:
##  $ ID            : chr  "GO:0051607" "GO:0009615" "GO:0045071" "GO:0045069" ...
##  $ Description   : chr  "defense response to virus" "response to virus" "negative regulation of viral genome replication" "regulation of viral genome replication" ...
##  $ GeneRatio     : chr  "5/6" "5/6" "3/6" "3/6" ...
##  $ BgRatio       : chr  "329/28928" "412/28928" "65/28928" "100/28928" ...
##  $ RichFactor    : num  0.0152 0.0121 0.0462 0.03 0.0261 ...
##  $ FoldEnrichment: num  73.3 58.5 222.5 144.6 125.8 ...
##  $ zScore        : num  19 16.9 25.8 20.7 19.3 ...
##  $ pvalue        : num  1.10e-09 3.39e-09 2.16e-07 7.96e-07 1.21e-06 ...
##  $ p.adjust      : num  1.37e-07 2.12e-07 8.98e-06 2.49e-05 3.03e-05 ...
##  $ qvalue        : num  3.35e-08 5.18e-08 2.19e-06 6.07e-06 7.41e-06 ...
##  $ geneID        : chr  "Oasl1/Rsad2/Ifit1/Oasl2/Rtp4" "Oasl1/Rsad2/Ifit1/Oasl2/Rtp4" "Oasl1/Rsad2/Oasl2" "Oasl1/Rsad2/Oasl2" ...
##  $ Count         : int  5 5 3 3 3 3 3 2 3 3 ...
## #...Citation
## G Yu. Thirteen years of clusterProfiler. The Innovation. 2024, 5(6):100722
 dotplot(pos_go) +
    ggtitle("") +
    theme_classic() + 
    theme(
        plot.title = element_text(hjust = 0.5),
        legend.position = "left",
        axis.text.y = element_text(hjust = 0, size = 10)) +
    scale_y_discrete(position = "right", 
                     labels = function(x) str_wrap(x, width = 25))  # Wrap y-axis labels to 2
## Scale for y is already present.
## Adding another scale for y, which will replace the existing scale.